# Figures 1, 2 and 3 for the paper "The Company You Keep: International Socialization and the Diffusion of Human Rights Norms"

# Please contact Brian Greenhill (bdgreen@u.washington.edu) with any questions regarding this code.

### Figure 1:
setwd("/Users/Brian/Brian/UW/Research Projects/IGOs/COW IGO data")
igodata<-read.csv("IGO_stateunit_v2.1.csv") # This file needs to be downloaded from the Correlates of War website (http://www.correlatesofwar.org).

# first get rid of the last column:
igodata<-igodata[,-500]

# now replace the 2,3,-9,-1 designations with zeros:
for (i in 5:499) {
	xvector<-igodata[,i]
	xvector[xvector==2|xvector==3|xvector==-9|xvector==-1]<-0
	igodata[,i]<-xvector
	print(i)
	}

# now create a vector of all the years in the dataset

yearlist<-sort(unique(igodata$year))

totaligos<-rep(NA,length(yearlist))

for (i in 1:length(yearlist)){
	
	# pick a year
	xyear<-yearlist[i]
	
	# pick out all the rows from the igodata matrix corresponding to that year:
	xmatrix<-igodata[igodata$year==xyear,]
	
	# now count up the number of igos present that year:
	coltotal<-rep(NA,499)
	for (a in 5:499){
		coltotal[a]<-sum(xmatrix[,a],na.rm=T)
		} # close a loop
		
		# now convert the coltotals into simply zeros or ones, to indicate presence or absence of that IGO.
		coltotal[coltotal>0]<-1
		
		# now sum up the total igos for that year:
		totaligos[i]<-sum(coltotal,na.rm=T)
		
		print(xyear)
		} # close the i loop

pdf(file="Figure 1.pdf", width=7, height=7)
plot(yearlist,totaligos,type="n",xlab="Year",xlim=c(1800,2000),ylab="Number of IGOs", las=1, bty="n")
lines(yearlist,totaligos,lwd=3)
dev.off()

### Figure 2:

setwd("/Users/Brian/Brian/UW/Research Projects/IGOs/Paper Drafts/ISQ Acceptance/ISQ Replication Files")
HR<-read.csv("Greenhill ISQ Replication dataset.csv")
attach(HR)

polity<-polity2
fdi<-fdi.in
log.gdp<-log(gdp)


# First create the one-year lagged physint variable:
physint1<-rep(NA,nrow(HR))

for (i in 2:nrow(HR)){
	if (bri.code[i]==bri.code[i-1]){
		physint1[i]<-physint[i-1]
	} # close the if condition
	} # close the i loop

# now create the 3-year lagged variables:
log.gdp3<-rep(NA,nrow(HR))
exports3<-rep(NA,nrow(HR))
fdi3<-rep(NA,nrow(HR))
fdistock3<-rep(NA,nrow(HR))
polity3<-rep(NA,nrow(HR))
bil.trade3<-rep(NA,nrow(HR))
language3<-rep(NA,nrow(HR))
colonies3<-rep(NA,nrow(HR))
neighbour3<-rep(NA,nrow(HR))
igo3<-rep(NA,nrow(HR))
ingo3<-rep(NA,nrow(HR))
hardpta3<-rep(NA,nrow(HR))
softpta3<-rep(NA,nrow(HR))
fdi.out3<-rep(NA,nrow(HR))
trade3<-rep(NA,nrow(HR))
density3<-rep(NA,nrow(HR))
durable3<-rep(NA,nrow(HR))
civwar3<-rep(NA,nrow(HR))
intwar3<-rep(NA,nrow(HR))
IGOc3<-rep(NA,nrow(HR))

for (i in 4:nrow(HR)){
	if (bri.code[i]==bri.code[i-3]){
		log.gdp3[i]<-log.gdp[i-3]
		exports3[i]<-exports[i-3]
		fdi3[i]<-fdi[i-3]
		fdistock3[i]<-fdistock[i-3]
		polity3[i]<-polity[i-3]
		bil.trade3[i]<-bil.trade[i-3]
		language3[i]<-language[i-3]
		colonies3[i]<-colonies[i-3]
		neighbour3[i]<-neighbour[i-3]
		igo3[i]<-igo[i-3]
		ingo3[i]<-ingo[i-3]
		hardpta3[i]<-hardpta[i-3]
		softpta3[i]<-softpta[i-3]
		fdi.out3[i]<-fdi.out[i-3]
		trade3[i]<-trade[i-3]
		density3[i]<-density[i-3]
		durable3[i]<-durable[i-3]
		civwar3[i]<-civwar[i-3]
		intwar3[i]<-intwar[i-3]
		IGOc3[i]<-IGOc[i-3]
		
		} # close the if condition
	} # close the i loop

library(MASS)
model1<-polr(as.factor(physint)~physint1+polity3+trade3+fdi3+log.gdp3+civwar3+intwar3+durable3+density3+hardpta3+softpta3+neighbour3+language3+colonies3+IGOc3,Hess=T,method="probit")

# now write a function to conduct the simulation:
simulate.IGOc<-function(vcmatrix=solve(model1$Hessian[1:15,1:15]), betas=coef(model1), cutpoints= model1$zeta, IGOc.range=seq(from=3.7,to=6.3,length.out=9), Xphysint1=median(physint1,na.rm=T), Xpolity3=median(polity3,na.rm=T), Xtrade3=median(trade3,na.rm=T), Xfdi3=median(fdi3,na.rm=T), Xlog.gdp3=median(log.gdp3,na.rm=T), Xcivwar3=median(civwar3,na.rm=T), Xintwar3=median(intwar3,na.rm=T), Xdurable3=median(durable3,na.rm=T), Xdensity3=median(density3,na.rm=T), Xhardpta3=median(hardpta3,na.rm=T), Xsoftpta3=median(softpta3,na.rm=T), Xneighbour3=median(neighbour3,na.rm=T), Xlanguage3=median(language3,na.rm=T), Xcolonies3=median(colonies3,na.rm=T), heading=""){

	mu<-matrix(NA,nrow=1000,ncol=9)
	mu.means<-rep(NA,9)
	mu.lowerCI<-rep(NA,9)
	mu.upperCI<-rep(NA,9)
	pi.hat<-matrix(NA,nrow=9,ncol=9)
	pi.lowerCI<-matrix(NA,nrow=9,ncol=9)
	pi.upperCI<-matrix(NA,nrow=9,ncol=9)

	for (i in 1:9) { # go through the range of IGO.c values
	
		# define xc:
		xc<-c(Xphysint1, Xpolity3, Xtrade3, Xfdi3, Xlog.gdp3, Xcivwar3, Xintwar3, Xdurable3, Xdensity3, Xhardpta3, Xsoftpta3, Xneighbour3, Xlanguage3, Xcolonies3, IGOc.range[i])
	
		# now calculate the expected value for a given draw:
		beta.draw<-mvrnorm(1000,betas,vcmatrix)
		mu[,i]<-beta.draw %*% xc
		mu.means[i]<-mean(mu[,i])
		mu.lowerCI[i]<-(sort(mu[,i]))[25]
		mu.upperCI[i]<-(sort(mu[,i]))[975]

		# now calculate the probability of each category
		pi.hat[1,i]<-pnorm(cutpoints[1],mu.means[i])
		pi.hat[2,i]<-pnorm(cutpoints[2],mu.means[i])-pnorm(cutpoints[1],mu.means[i])
		pi.hat[3,i]<-pnorm(cutpoints[3],mu.means[i])-pnorm(cutpoints[2],mu.means[i])
		pi.hat[4,i]<-pnorm(cutpoints[4],mu.means[i])-pnorm(cutpoints[3],mu.means[i])
		pi.hat[5,i]<-pnorm(cutpoints[5],mu.means[i])-pnorm(cutpoints[4],mu.means[i])
		pi.hat[6,i]<-pnorm(cutpoints[6],mu.means[i])-pnorm(cutpoints[5],mu.means[i])
		pi.hat[7,i]<-pnorm(cutpoints[7],mu.means[i])-pnorm(cutpoints[6],mu.means[i])
		pi.hat[8,i]<-pnorm(cutpoints[8],mu.means[i])-pnorm(cutpoints[7],mu.means[i])
		pi.hat[9,i]<-1-pnorm(cutpoints[8],mu.means[i])
		
		} # close i loop
		
	# now convert the pi.hat matrix into a cumulative measure
	pi.hat2<-matrix(NA,nrow=9,ncol=9)
	pi.hat2[1,]<-pi.hat[1,]
	for (i in 2:9) {
		pi.hat2[i,]<-pi.hat[i,]+pi.hat2[i-1,]
		}

		plot(IGOc.range,0:8,type="n",ylim=c(0,1),xlab="IGO Context",ylab="Cumulative Probability of Response",main=heading)

		polygon(x=c(IGOc.range,rev(IGOc.range)),y=c(rep(0,9),rev(pi.hat2[1,])),col="gray10")

		col.range<-c("gray20","gray30","gray40","gray50","gray65","gray80","gray90","gray95")

		for (i in 1:8) {
			polygon(x=c(IGOc.range,rev(IGOc.range)),y=c(pi.hat2[i,],rev(pi.hat2[i+1,])),col=col.range[i])
			}
			
		# add text markers
		xrange<-seq(from=IGOc.range[2], to=IGOc.range[8], length.out=9)
		for (i in 1:8) {
			ylevel<-mean(c(pi.hat2[i,i+1],pi.hat2[i+1,i+1]))
			if (ylevel>0.98) ylevel<-0.98
			if (ylevel<0.02) ylevel<-0.02
			text(xrange[i+1],ylevel,i,cex=1.3)
			}
		# do the zero label:
		ylevel<-mean(c(0,pi.hat2[1,1]))
		if (ylevel<0.02) ylevel<-0.02
		text(xrange[1],ylevel,0,cex=1.3)

		for (i in 1:9) {
			lines(IGOc.range,pi.hat2[i,],lwd=2)
			}	

	} # close simulate.IGOc function

vcov.stata<-read.csv("/Users/Brian/Brian/UW/Research Projects/IGOs/Other data/varcovar.csv")
vcov.stata<-vcov.stata[1:15,2:16]

# Now run a few simulations:
pdf(file="Figure 2.pdf", width=8, height=4)
par(mfrow=c(1,3))
simulate.IGOc(heading="Median State \n (GDP= $1,700; Polity= 0)", vcmatrix=vcov.stata)
simulate.IGOc(Xpolity3=quantile(polity3,0.1,na.rm=T),Xlog.gdp3=quantile(log.gdp3,0.1,na.rm=T), heading="Poor Authoritarian State \n (GDP= $240; Polity= -8)", vcmatrix=vcov.stata)
simulate.IGOc(Xpolity3=quantile(polity3,0.9,na.rm=T),Xlog.gdp3=quantile(log.gdp3,0.9,na.rm=T), heading="Rich Democratic State \n (GDP= $18,000; Polity= +10)", vcmatrix=vcov.stata)
dev.off()




### Figure 3:

x1<-1:5
x2<-x1+0.15
# IGOc:
b1<-c(0.256,0.227,0.375,0.231,0.207)
l1<-c(0.069,0.028,0.169,0.063,0.012)
u1<-c(0.444,0.426,0.582,0.399,0.402)
# IGOcxhr:
b2<-c(0.254,0.227,0.367,0.227,0.204)
l2<-c(0.069,0.030,0.163,0.061,0.011)
u2<-c(0.440,0.423,0.571,0.394,0.398)

pdf(file="Figure 3.pdf", width=7, height=7)
par(mfrow=c(1,1))
plot(x1,b1,type="n",xlim=c(1,5.2),ylim=c(-0.2,0.8),xlab="Time Lag (years)",ylab="Coefficient of IGO Context", las=1, bty="n")
abline(h=0,lty=2)

lines(x1,b1,lwd=3,col="grey50")
points(x1,b1,pch=21,cex=1.4,col="grey50",bg="grey50")

for (i in 1:5) {
	lines(x=c(x1[i],x1[i]),y=c(l1[i],u1[i]),col="grey50")
	lines(x=c(x1[i]-0.05,x1[i]+0.05),y=c(l1[i],l1[i]),col="grey50")
	lines(x=c(x1[i]-0.05,x1[i]+0.05),y=c(u1[i],u1[i]),col="grey50")
	}
	
lines(x2,b2,lwd=3)
points(x2,b2,pch=22,cex=1.4,bg=1)

for (i in 1:5) {
	lines(x=c(x2[i],x2[i]),y=c(l2[i],u2[i]))
	lines(x=c(x2[i]-0.05,x2[i]+0.05),y=c(l2[i],l2[i]))
	lines(x=c(x2[i]-0.05,x2[i]+0.05),y=c(u2[i],u2[i]))
	}
	
dev.off()


